home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
rwvector.lha
/
RWVector2.1
/
src
/
xludecmp.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1989-08-18
|
3KB
|
121 lines
/*
* <T> Precision LU Decomposition
*
* Copyright (C) 1988, 1989.
*
* Dr. Thomas Keffer
* Rogue Wave Associates
* P.O. Box 85341
* Seattle WA 98145-1341
*
* Permission to use, copy, modify, and distribute this
* software and its documentation for any purpose and
* without fee is hereby granted, provided that the
* above copyright notice appear in all copies and that
* both that copyright notice and this permission notice
* appear in supporting documentation.
*
* This software is provided "as is" without any
* expressed or implied warranty.
*
*
* @(#)xludecmp.cc 2.1 8/18/89
*/
/*
* These routines call the appropriate LINPACK routines
*/
#define NO_VECTOR_MATHFUN
#include "rw/<A>LUDecomp.h"
#include "rw/linpack.h"
#include <stdio.h>
// Construct an LU decomposition from m
<A>LUDecomp::<A>LUDecomp(const <A>GEMatrix& m)
: (m.copy())
{
assertSquare();
order = m.rows();
ipvt = new fortran_int[order];
<A>gefa_(data(), &order, &order, ipvt, &info);
}
// Get the determinant from an existing LU decomposition
<T>
determinant(const <A>LUDecomp& d)
{
d.assertDefined();
d.assertPivots();
<T> determine[2];
<T>* temp= new <T>[d.order];
fortran_int job = 10; // Indicate only determinant to be computed
<A>gedi_(d.data(), &d.order, &d.order, d.ipvt, determine,
temp, &job);
delete temp;
return determine[1] * pow(10.0, determine[2]);
}
// Calculate the inverse of an existing LU decomposition
<A>GEMatrix
inverse(const <A>LUDecomp& d)
{
d.assertDefined();
d.assertPivots();
<T> determine[2];
<T>* temp= new <T>[d.order];
fortran_int job = 1; // Indicate only inverse to be computed
<A>GEMatrix the_inverse = d.copy();
<A>gedi_(the_inverse.data(), &d.order, &d.order, d.ipvt, determine,
temp, &job);
delete temp;
return the_inverse;
}
<T>Vec
solve(const <A>LUDecomp& d, const <T>Vec& rhs)
{
d.assertDefined();
d.assertPivots();
<T>Vec temp = rhs.copy();
fortran_int job = 0;
<A>gesl_(d.data(), &d.order, &d.order, d.ipvt, temp.data(), &job);
return temp;
}
/************ Utility routines ******************/
void
<A>LUDecomp::assertDefined()
{
// Check to make sure this isn't a null matrix
if(!order){
char msg[40];
sprintf(msg, "Undefined LU Decomposition");
RWnote("<A>LUDecomp::assertDefined()", msg);
RWerror(DEFAULT);
}
}
void
<A>LUDecomp::assertPivots()
{
// Check to make sure none of the pivots are zero
if(info){
char msg[120];
sprintf(msg, "Pivot number %d is zero", info);
RWnote("<A>LUDecomp::assertPivot()", msg);
RWerror(DEFAULT);
}
}